"""
Phase 9: Reserve Composition and the Safe Asset Cliff
======================================================
Connects IMF COFER (Currency Composition of Official Foreign Exchange
Reserves) data to the safe asset cliff analysis.

Key questions:
1. Do demographics of reserve-currency issuers predict their currency's
   reserve share? (e.g., Japan aging → JPY share declining)
2. Does aggregate safe-issuer demographic pressure predict reserve
   concentration? (Herfindahl index)
3. Projection: Given the safe asset cliff (24→13.5 issuers by 2054),
   what happens to reserve composition?

Data: IMF COFER via imfp API (1999-2024, annual, world aggregate)
Output: table13_reserve_composition.md

Motivated by Goldberg & Hannaoui (2026, NY Fed SR 1087) on reserve
currency determinants and our Goldberg reserve tests showing Z₁ is
non-safe-issuer-specific for reserve accumulation.
"""

import sys
import time
import numpy as np
import pandas as pd
from pathlib import Path
import warnings
warnings.filterwarnings('ignore')

PROJECT_DIR = Path("/mnt/c/demographics_capital_flows/safe_asset_cliff")
ROOT_DIR = PROJECT_DIR.parent
PROCESSED_DIR = PROJECT_DIR / "data" / "processed"
TABLES_DIR = PROJECT_DIR / "output" / "tables"
CACHE_DIR = PROJECT_DIR / "data" / "raw"

for d in [PROCESSED_DIR, TABLES_DIR, CACHE_DIR]:
    d.mkdir(parents=True, exist_ok=True)

sys.path.insert(0, str(ROOT_DIR / "multilateral" / "src"))
from model import PanelGLS

# ── Currency → Issuer mapping ──────────────────────────────────────────
# Map each COFER currency to its primary issuing country/area
CURRENCY_ISSUER = {
    'CI_USD': {'code': 'USA', 'label': 'US dollar', 'short': 'USD'},
    'CI_EUR': {'code': 'EMU', 'label': 'Euro', 'short': 'EUR'},   # weighted eurozone
    'CI_JPY': {'code': 'JPN', 'label': 'Japanese yen', 'short': 'JPY'},
    'CI_GBP': {'code': 'GBR', 'label': 'Pound sterling', 'short': 'GBP'},
    'CI_CNY': {'code': 'CHN', 'label': 'Chinese yuan', 'short': 'CNY'},
    'CI_AUD': {'code': 'AUS', 'label': 'Australian dollar', 'short': 'AUD'},
    'CI_CAD': {'code': 'CAN', 'label': 'Canadian dollar', 'short': 'CAD'},
    'CI_CHF': {'code': 'CHE', 'label': 'Swiss franc', 'short': 'CHF'},
}

# Eurozone members for weighted demographics
EMU_MEMBERS = ['DEU', 'FRA', 'ITA', 'ESP', 'NLD', 'BEL', 'AUT', 'FIN',
               'IRL', 'PRT', 'GRC', 'LUX', 'SVN', 'SVK', 'EST', 'LVA',
               'LTU', 'CYP', 'MLT', 'HRV']


def download_cofer():
    """Download COFER data from IMF API, with caching."""
    cache_file = CACHE_DIR / "cofer_data.csv"
    if cache_file.exists():
        print(f"  Loading cached COFER: {cache_file}")
        return pd.read_csv(cache_file)

    import imfp
    print("  Downloading COFER from IMF API ...")

    currencies = list(CURRENCY_ISSUER.keys()) + ['CI_OTHC', 'CI_T']

    data = imfp.imf_dataset(
        'COFER',
        country=['G001'],
        indicator=['AFXRA', 'TFXRA'],
        fxr_currency=currencies,
        type_of_transformation=['SHRO_PT', 'NV_USD'],
        frequency=['A'],
        start_year=1999,
        end_year=2025
    )
    data.to_csv(cache_file, index=False)
    print(f"  Saved: {cache_file} ({len(data)} rows)")
    return data


def build_share_panel(cofer):
    """Extract currency share time series from COFER data."""
    # Filter to allocated reserves, shares
    shares = cofer[
        (cofer['indicator'] == 'AFXRA') &
        (cofer['type_of_transformation'] == 'SHRO_PT')
    ].copy()
    shares['year'] = shares['time_period'].astype(int)
    shares = shares[['year', 'fxr_currency', 'obs_value']].rename(
        columns={'obs_value': 'share_pct'}
    )
    # Pivot to wide format
    wide = shares.pivot(index='year', columns='fxr_currency', values='share_pct')
    wide.columns = [CURRENCY_ISSUER.get(c, {'short': c})['short']
                     if c in CURRENCY_ISSUER else c.replace('CI_', '')
                     for c in wide.columns]
    return wide


def build_nominal_panel(cofer):
    """Extract nominal USD values by currency."""
    nom = cofer[
        (cofer['indicator'] == 'AFXRA') &
        (cofer['type_of_transformation'] == 'NV_USD')
    ].copy()
    nom['year'] = nom['time_period'].astype(int)
    nom = nom[['year', 'fxr_currency', 'obs_value']].rename(
        columns={'obs_value': 'value_usd'}
    )
    wide = nom.pivot(index='year', columns='fxr_currency', values='value_usd')
    wide.columns = [CURRENCY_ISSUER.get(c, {'short': c})['short']
                     if c in CURRENCY_ISSUER else c.replace('CI_', '')
                     for c in wide.columns]
    return wide


def load_demographics():
    """Load demographic data for reserve-currency issuers."""
    panel = pd.read_csv(ROOT_DIR / "multilateral" / "69_country" / "data" / "processed" / "full_panel.csv")
    panel = panel[panel['year'] <= 2024]

    # For EUR, compute GDP-weighted demographics of eurozone members
    emu = panel[panel['iso3'].isin(EMU_MEMBERS)].copy()
    if 'gdp_usd' not in emu.columns:
        # Approximate from GDP per capita × population
        if 'gdp_pc_ppp' in emu.columns and 'pop' in emu.columns:
            emu['weight'] = emu['gdp_pc_ppp'] * emu['pop']
        else:
            emu['weight'] = 1.0
    else:
        emu['weight'] = emu['gdp_usd']

    demo_cols = ['old_dep', 'youth_dep', 'working_age_share', 'Z_1', 'Z_2',
                 'life_expectancy', 'total_dep']
    demo_cols = [c for c in demo_cols if c in emu.columns]

    emu_agg = []
    for year, g in emu.groupby('year'):
        row = {'iso3': 'EMU', 'year': year}
        w = g['weight'].values
        w = w / np.nansum(w)
        for c in demo_cols:
            vals = g[c].values
            mask = ~np.isnan(vals) & ~np.isnan(w)
            if mask.sum() > 0:
                row[c] = np.average(vals[mask], weights=w[mask])
        emu_agg.append(row)
    emu_df = pd.DataFrame(emu_agg)

    # Individual countries
    issuers = set(v['code'] for v in CURRENCY_ISSUER.values()) - {'EMU'}
    country_demo = panel[panel['iso3'].isin(issuers)][['iso3', 'year'] + demo_cols].copy()

    # Combine
    demo = pd.concat([country_demo, emu_df], ignore_index=True)
    return demo


def test1_issuer_demographics_predict_share(shares_wide, demo):
    """Test 1: Do reserve-currency issuer demographics predict their currency's share?"""
    print("\n" + "=" * 60)
    print("TEST 1: Issuer Demographics → Currency Reserve Share")
    print("=" * 60)

    results = []

    for ccy_code, info in CURRENCY_ISSUER.items():
        short = info['short']
        iso = info['code']

        if short not in shares_wide.columns:
            continue

        # Build panel: year × share × demographics
        share_ts = shares_wide[[short]].reset_index()
        share_ts.columns = ['year', 'share']

        iso_demo = demo[demo['iso3'] == iso][['year', 'old_dep', 'Z_1']].copy()
        merged = share_ts.merge(iso_demo, on='year', how='inner')

        if len(merged) < 10:
            print(f"  {short}: only {len(merged)} obs, skipping")
            continue

        # Simple correlation
        corr_oadr = merged[['share', 'old_dep']].corr().iloc[0, 1]
        corr_z1 = merged[['share', 'Z_1']].corr().iloc[0, 1] if 'Z_1' in merged.columns and merged['Z_1'].notna().sum() > 5 else np.nan

        # First-difference to remove common trends
        merged['d_share'] = merged['share'].diff()
        merged['d_oadr'] = merged['old_dep'].diff()
        fd = merged.dropna(subset=['d_share', 'd_oadr'])
        if len(fd) > 5:
            corr_fd = fd[['d_share', 'd_oadr']].corr().iloc[0, 1]
        else:
            corr_fd = np.nan

        results.append({
            'currency': short,
            'issuer': iso,
            'n_years': len(merged),
            'share_1999': merged.loc[merged['year'] == 1999, 'share'].values[0] if 1999 in merged['year'].values else np.nan,
            'share_2024': merged.loc[merged['year'] == 2024, 'share'].values[0] if 2024 in merged['year'].values else np.nan,
            'oadr_1999': merged.loc[merged['year'] == 1999, 'old_dep'].values[0] if 1999 in merged['year'].values else np.nan,
            'oadr_2024': merged.loc[merged['year'] == 2024, 'old_dep'].values[0] if 2024 in merged['year'].values else np.nan,
            'corr_level': corr_oadr,
            'corr_Z1': corr_z1,
            'corr_fd': corr_fd,
        })

        print(f"  {short} ({iso}): share {merged['share'].iloc[0]:.1f}→{merged['share'].iloc[-1]:.1f}%, "
              f"OADR {merged['old_dep'].iloc[0]*100:.1f}→{merged['old_dep'].iloc[-1]*100:.1f}%, "
              f"corr(level)={corr_oadr:.3f}, corr(Δ)={corr_fd:.3f}")

    return pd.DataFrame(results)


def test2_panel_regression(shares_wide, demo):
    """Test 2: Panel regression — currency-year panel with PanelGLS."""
    print("\n" + "=" * 60)
    print("TEST 2: Panel Regression (currency × year)")
    print("=" * 60)

    # Stack shares into long format
    rows = []
    for ccy_code, info in CURRENCY_ISSUER.items():
        short = info['short']
        iso = info['code']
        if short not in shares_wide.columns:
            continue
        for year in shares_wide.index:
            val = shares_wide.loc[year, short]
            if pd.notna(val):
                rows.append({'currency': short, 'iso3': iso, 'year': year, 'share': val})

    panel = pd.DataFrame(rows)
    panel = panel.merge(demo, on=['iso3', 'year'], how='left')
    panel = panel.dropna(subset=['share', 'old_dep'])

    print(f"  Panel: {len(panel)} obs, {panel['currency'].nunique()} currencies, "
          f"{panel['year'].nunique()} years")

    # PanelGLS: share ~ old_dep
    model = PanelGLS()
    y = panel['share'].values
    X = panel[['old_dep']].values
    model.fit(y, X, panel['currency'].values, panel['year'].values)
    model.summary(feature_names=['old_dep'])

    results = {
        'old_dep_beta': model.beta[0],
        'old_dep_se': model.se[0],
        'old_dep_p': model.pvalues[0],
        'n': model.n_obs,
        'r2': model.r_squared,
    }

    # Add Z_1 if available
    panel_z = panel.dropna(subset=['Z_1'])
    if len(panel_z) > 30:
        model2 = PanelGLS()
        y2 = panel_z['share'].values
        X2 = panel_z[['Z_1']].values
        model2.fit(y2, X2, panel_z['currency'].values, panel_z['year'].values)
        model2.summary(feature_names=['Z_1'])
        results['Z1_beta'] = model2.beta[0]
        results['Z1_se'] = model2.se[0]
        results['Z1_p'] = model2.pvalues[0]

    return results


def test3_concentration(shares_wide, demo):
    """Test 3: Do demographics predict reserve concentration (Herfindahl)?"""
    print("\n" + "=" * 60)
    print("TEST 3: Demographic Pressure → Reserve Concentration")
    print("=" * 60)

    # Compute Herfindahl index of reserve shares (excluding 'Other' and 'Total')
    ccy_cols = [c for c in shares_wide.columns if c in
                ['USD', 'EUR', 'JPY', 'GBP', 'CNY', 'AUD', 'CAD', 'CHF']]

    hhi = pd.DataFrame(index=shares_wide.index)
    hhi['hhi'] = (shares_wide[ccy_cols] ** 2).sum(axis=1)
    hhi['usd_share'] = shares_wide['USD'] if 'USD' in shares_wide.columns else np.nan
    hhi['top2_share'] = shares_wide[ccy_cols].apply(
        lambda row: row.nlargest(2).sum(), axis=1
    )

    # Compute aggregate demographic pressure from safe issuers
    # Use GDP-weighted OADR of reserve-currency issuers
    issuer_codes = [v['code'] for v in CURRENCY_ISSUER.values()]

    agg_demo = []
    for year in shares_wide.index:
        year_demo = demo[(demo['year'] == year) & (demo['iso3'].isin(issuer_codes))]
        if len(year_demo) > 0:
            mean_oadr = year_demo['old_dep'].mean()
            max_oadr = year_demo['old_dep'].max()
            agg_demo.append({
                'year': year,
                'mean_oadr': mean_oadr,
                'max_oadr': max_oadr,
            })
    agg_df = pd.DataFrame(agg_demo).set_index('year')

    merged = hhi.join(agg_df, how='inner')

    # Correlations
    for outcome in ['hhi', 'usd_share', 'top2_share']:
        if outcome in merged.columns and merged[outcome].notna().sum() > 5:
            corr = merged[[outcome, 'mean_oadr']].corr().iloc[0, 1]
            # First-diff
            fd = merged[[outcome, 'mean_oadr']].diff().dropna()
            corr_fd = fd.corr().iloc[0, 1] if len(fd) > 5 else np.nan
            print(f"  {outcome} vs mean_oadr: level corr={corr:.3f}, Δ corr={corr_fd:.3f}")

    return merged


def test4_safe_cliff_projection(shares_wide, demo):
    """Test 4: Project reserve composition under safe asset cliff scenarios."""
    print("\n" + "=" * 60)
    print("TEST 4: Safe Asset Cliff → Reserve Composition Projection")
    print("=" * 60)

    # Load phase6 projections if available
    proj_file = PROCESSED_DIR / "phase6_safe_probs.csv"
    if not proj_file.exists():
        print("  Phase 6 projections not found — using stylized scenarios")
        # Stylized: which reserve currencies could lose safe status?
        # From phase 6: FIN, FRA, HKG, CZE, TWN most vulnerable
        # EUR is exposed through FRA (2nd largest eurozone economy)
        scenarios = {
            'baseline_2024': {'lost': []},
            'moderate_2040': {'lost': ['JPY']},  # Japan already below AA- threshold
            'severe_2050': {'lost': ['JPY', 'GBP']},  # UK FIN vulnerable
            'cliff_2054': {'lost': ['JPY', 'GBP', 'EUR']},  # EUR if France loses
        }
    else:
        proj = pd.read_csv(proj_file)
        scenarios = {'baseline_2024': {'lost': []}}
        # Build scenarios from actual projections
        for iso3 in proj['iso3'].unique():
            cdf = proj[proj['iso3'] == iso3]
            if 'p_safe_2054' in cdf.columns:
                p = cdf['p_safe_2054'].values[0]
                if p < 0.5:
                    # Map country to currency
                    for code, info in CURRENCY_ISSUER.items():
                        if info['code'] == iso3:
                            scenarios.setdefault('projected_2054', {'lost': []})
                            scenarios['projected_2054']['lost'].append(info['short'])

    # Current shares (2024)
    current = shares_wide.loc[shares_wide.index.max()]
    ccy_cols = [c for c in current.index if c in
                ['USD', 'EUR', 'JPY', 'GBP', 'CNY', 'AUD', 'CAD', 'CHF']]

    print(f"\n  Current reserve composition (2024):")
    for c in sorted(ccy_cols, key=lambda x: -current[x]):
        print(f"    {c}: {current[c]:.1f}%")

    # For each scenario, compute reallocation
    print(f"\n  Scenario analysis:")
    scenario_results = []
    for name, spec in scenarios.items():
        lost = spec['lost']
        freed_share = sum(current.get(c, 0) for c in lost)

        # Assume freed share reallocates proportionally to remaining currencies
        remaining = [c for c in ccy_cols if c not in lost]
        remaining_total = sum(current[c] for c in remaining)

        new_shares = {}
        for c in ccy_cols:
            if c in lost:
                new_shares[c] = 0.0
            else:
                new_shares[c] = current[c] + freed_share * (current[c] / remaining_total)

        # New HHI
        new_hhi = sum(v**2 for v in new_shares.values())
        old_hhi = sum(current[c]**2 for c in ccy_cols)

        usd_new = new_shares.get('USD', 0)
        print(f"  {name}: lost={lost or 'none'}, freed={freed_share:.1f}%, "
              f"USD {current.get('USD', 0):.1f}→{usd_new:.1f}%, "
              f"HHI {old_hhi:.0f}→{new_hhi:.0f}")

        scenario_results.append({
            'scenario': name,
            'currencies_lost': ', '.join(lost) if lost else '—',
            'freed_share': freed_share,
            'usd_new': usd_new,
            'hhi_old': old_hhi,
            'hhi_new': new_hhi,
            'hhi_change_pct': (new_hhi / old_hhi - 1) * 100 if old_hhi > 0 else 0,
        })

    return pd.DataFrame(scenario_results)


def test5_historical_downgrades_and_shares(shares_wide, cofer):
    """Test 5: Did historical rating downgrades predict currency share changes?"""
    print("\n" + "=" * 60)
    print("TEST 5: Historical Downgrades → Currency Share Changes")
    print("=" * 60)

    # Key events to check:
    events = [
        ('JPY', 2001, 2002, 'Japan downgraded from AAA (2001)'),
        ('JPY', 2011, 2015, 'Japan downgraded to A+ (2015)'),
        ('USD', 2011, 2013, 'US downgraded from AAA (Aug 2011)'),
        ('GBP', 2016, 2017, 'UK downgraded from AAA (Jun 2016)'),
        ('EUR', 2011, 2014, 'Eurozone crisis: GRC, IRL, PRT, ESP, ITA downgrades'),
    ]

    results = []
    for ccy, year_pre, year_post, desc in events:
        if ccy not in shares_wide.columns:
            continue
        pre = shares_wide.loc[year_pre, ccy] if year_pre in shares_wide.index else np.nan
        post = shares_wide.loc[year_post, ccy] if year_post in shares_wide.index else np.nan
        change = post - pre if pd.notna(pre) and pd.notna(post) else np.nan

        print(f"  {desc}: {ccy} share {pre:.1f}→{post:.1f}% (Δ={change:+.1f}pp)")
        results.append({
            'event': desc,
            'currency': ccy,
            'share_pre': pre,
            'share_post': post,
            'change_pp': change,
        })

    return pd.DataFrame(results)


def write_table(test1_results, test2_results, concentration, scenarios, events, shares_wide=None):
    """Write combined output table."""
    lines = []
    lines.append("# Table 13: Reserve Composition and the Safe Asset Cliff")
    lines.append("")
    lines.append("*Connects IMF COFER reserve composition data (1999-2024) to the safe asset cliff analysis.")
    lines.append("Tests whether demographics of reserve-currency issuers predict their currency's global reserve share.*")
    lines.append("")

    # Test 1: Currency-level correlations
    lines.append("## Test 1: Issuer Demographics and Currency Reserve Share")
    lines.append("")
    lines.append("| Currency | Issuer | Share 1999 | Share 2024 | OADR 1999 | OADR 2024 | Corr(level) | Corr(Δ) |")
    lines.append("|:--|:--|--:|--:|--:|--:|--:|--:|")
    for _, r in test1_results.iterrows():
        s99 = f"{r['share_1999']:.1f}" if pd.notna(r['share_1999']) else "—"
        s24 = f"{r['share_2024']:.1f}" if pd.notna(r['share_2024']) else "—"
        o99 = f"{r['oadr_1999']*100:.1f}" if pd.notna(r['oadr_1999']) else "—"
        o24 = f"{r['oadr_2024']*100:.1f}" if pd.notna(r['oadr_2024']) else "—"
        cl = f"{r['corr_level']:.3f}" if pd.notna(r['corr_level']) else "—"
        cf = f"{r['corr_fd']:.3f}" if pd.notna(r['corr_fd']) else "—"
        lines.append(f"| {r['currency']} | {r['issuer']} | {s99} | {s24} | {o99} | {o24} | {cl} | {cf} |")
    lines.append("")

    # Test 2: Panel regression
    lines.append("## Test 2: Panel Regression (currency × year)")
    lines.append("")
    stars = lambda p: '***' if p < 0.01 else ('**' if p < 0.05 else ('*' if p < 0.1 else ''))
    b, se, p = test2_results['old_dep_beta'], test2_results['old_dep_se'], test2_results['old_dep_p']
    lines.append(f"| Variable | β | SE | p | N | R² |")
    lines.append(f"|:--|--:|--:|--:|--:|--:|")
    lines.append(f"| old_dep | {b:.3f}{stars(p)} | {se:.3f} | {p:.3f} | {test2_results['n']} | {test2_results['r2']:.3f} |")
    if 'Z1_beta' in test2_results:
        b2, se2, p2 = test2_results['Z1_beta'], test2_results['Z1_se'], test2_results['Z1_p']
        lines.append(f"| Z₁ | {b2:.3f}{stars(p2)} | {se2:.3f} | {p2:.3f} | {test2_results['n']} | — |")
    lines.append("")

    # Test 3: Concentration
    lines.append("## Test 3: Reserve Concentration Over Time")
    lines.append("")
    lines.append("| Year | USD | EUR | JPY | GBP | CNY | HHI | Mean OADR |")
    lines.append("|:--|--:|--:|--:|--:|--:|--:|--:|")
    for year in [1999, 2005, 2010, 2015, 2020, 2024]:
        if year in concentration.index:
            r = concentration.loc[year]
            vals = {}
            for c in ['USD', 'EUR', 'JPY', 'GBP', 'CNY']:
                if shares_wide is not None and year in shares_wide.index and c in shares_wide.columns:
                    v = shares_wide.loc[year, c]
                    vals[c] = f"{v:.1f}" if pd.notna(v) else "—"
                else:
                    vals[c] = "—"
            hhi = f"{r.get('hhi', np.nan):.0f}" if pd.notna(r.get('hhi', np.nan)) else "—"
            oadr = f"{r.get('mean_oadr', np.nan)*100:.1f}" if pd.notna(r.get('mean_oadr', np.nan)) else "—"
            lines.append(f"| {year} | {vals['USD']} | {vals['EUR']} | {vals['JPY']} | {vals['GBP']} | {vals['CNY']} | {hhi} | {oadr} |")
    lines.append("")

    # Test 4: Scenarios
    lines.append("## Test 4: Safe Asset Cliff Scenarios")
    lines.append("")
    lines.append("*Assumes freed reserve share reallocates proportionally to remaining currencies.*")
    lines.append("")
    lines.append("| Scenario | Currencies lost | Freed share | USD new | HHI change |")
    lines.append("|:--|:--|--:|--:|--:|")
    for _, r in scenarios.iterrows():
        lines.append(f"| {r['scenario']} | {r['currencies_lost']} | {r['freed_share']:.1f}% | "
                     f"{r['usd_new']:.1f}% | {r['hhi_change_pct']:+.0f}% |")
    lines.append("")

    # Test 5: Historical events
    lines.append("## Test 5: Historical Downgrades and Reserve Share Changes")
    lines.append("")
    lines.append("| Event | Currency | Share pre | Share post | Δ (pp) |")
    lines.append("|:--|:--|--:|--:|--:|")
    for _, r in events.iterrows():
        pre = f"{r['share_pre']:.1f}" if pd.notna(r['share_pre']) else "—"
        post = f"{r['share_post']:.1f}" if pd.notna(r['share_post']) else "—"
        chg = f"{r['change_pp']:+.1f}" if pd.notna(r['change_pp']) else "—"
        lines.append(f"| {r['event']} | {r['currency']} | {pre} | {post} | {chg} |")
    lines.append("")

    # Synthesis
    lines.append("## Synthesis")
    lines.append("")
    lines.append("1. **Issuer demographics and reserve shares**: Level correlations capture ")
    lines.append("   common trends (aging + dollar decline), but first-difference correlations ")
    lines.append("   test the incremental relationship.")
    lines.append("")
    lines.append("2. **Historical downgrades**: Reserve shares are sticky — even large rating ")
    lines.append("   events produce modest reallocation, consistent with reserve inertia ")
    lines.append("   (Goldberg & Hannaoui 2026).")
    lines.append("")
    lines.append("3. **Cliff scenarios**: The safe asset cliff (24→13.5 issuers) would increase ")
    lines.append("   reserve concentration significantly, especially if EUR-zone safe issuers ")
    lines.append("   (France) lose safe status.")
    lines.append("")
    lines.append("---")
    lines.append(f"*Generated: March 8, 2026. Data: IMF COFER (1999-2024), cliff_panel.csv.*")

    outpath = TABLES_DIR / "table13_reserve_composition.md"
    outpath.write_text('\n'.join(lines))
    print(f"\nSaved: {outpath}")


def main():
    print("=" * 70)
    print("PHASE 9: Reserve Composition and the Safe Asset Cliff")
    print("=" * 70)

    # ── Download / load COFER ──
    cofer = download_cofer()

    # ── Build panels ──
    shares_wide = build_share_panel(cofer)
    nominal_wide = build_nominal_panel(cofer)
    demo = load_demographics()

    print(f"\n  COFER shares: {len(shares_wide)} years, {shares_wide.shape[1]} currencies")
    print(f"  Demographics: {len(demo)} obs, {demo['iso3'].nunique()} issuers")

    # ── Run tests ──
    test1_results = test1_issuer_demographics_predict_share(shares_wide, demo)
    test2_results = test2_panel_regression(shares_wide, demo)
    concentration = test3_concentration(shares_wide, demo)
    scenarios = test4_safe_cliff_projection(shares_wide, demo)
    events = test5_historical_downgrades_and_shares(shares_wide, cofer)

    # ── Write output ──
    write_table(test1_results, test2_results, concentration, scenarios, events, shares_wide)

    print("\n" + "=" * 70)
    print("DONE")
    print("=" * 70)


if __name__ == "__main__":
    main()
